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Abstract 

Interactions between genes and gene products give rise to complex circuits that enable cells to pro- 
cess information and respond to external signals. Theoretical studies often describe these interactions 
using continuous, stochastic, or logical approaches. We propose a new modeling framework for gene 
regulatory networks, that combines the intuitive appeal of a qualitative description of gene states with a 
high flexibility in incorporating stochasticity in the duration of cellular processes. We apply our methods 
to the regulatory network of the segment polarity genes, thus gaining novel insights into the development 
of gene expression patterns. For example, we show that very short synthesis and decay times can perturb 
the wild type pattern. On the other hand, separation of timescales between pre- and posttranslational 
processes and a minimal prepattern ensure convergence to the wild type expression pattern regardless of 
fluctuations. 
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1 Introduction 

Understanding how genetic information is translated into proteins to produce various cell types remains a 
major challenge in contemporary biology |Wolpert et al, 1998) . Gene products often regulate the synthesis 
of mRNAs and proteins, forming complex networks of regulatory interactions. Concurrently with experi- 
mental progress in gene control networks [D avidson et al., 2002 |, several alternative modeling frameworks 
have been proposed. In the continuous-state approach, the concentrations of cellular components are as- 
sumed to be continuous functions of time, governed by differential equations with mass-action (or more 
general) kinetics |Reinitz and Sharp, 1995[|von Dassow et al., 2000| Gursky et al., 20 01 1. Stochastic mod- 



els address the deviations from population homogeneity by transforming reaction rates into probabilities and 
concentrations into numbers of molecules | Rao et al., 2002| . Finally, in the discrete approach, each compo- 
nent is assumed to have a small number of qualitative states, and the regulatory interactions are described by 
logical functions |Mendoza et al, 1999[|Sanchez and Thieffry, 2001||Yuh et al., 200Tl|Kauffman et al, 2003| 



Ghysen and Thomas, 2003j|Bodnar, 1997| [Albert and Othmer, 20031 . 



The kinetic details of protein-protein or protein-DNA interactions are rarely known, but there is in- 
creasing evidence that the input-output curves of regulatory relationships are strongly sigmoidal and can be 
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well approximated by step functions |Yuh et al., 200H|Thomas, 1973| . Moreover, both models and exper- 
iments suggest that regulatory networks are remarkably robust, that is, they maintain their function even 
when faced with fluctuations in components and reaction rates [ von D assow et al, 2000[|Alon et al, 1999| 
lEldar et al., 2002| Carlson and Doyle, 2002 Conant and Wagner, 2004J . These observations lend support to 



the assumption of discrete states for genetic network components and of combinatorial rules for the effects 
of transcription factors [Glass anS^auffman, 1973| de Jong et al., 2004] . The extreme of discretization. 



Boolean models, consider only two states (expressed or not), closely mimicking the inference methods used 
in genetics [Kauffman et al, 2003l|Thomas, 1973||Kauffman, 1993| . It is straightforward to study the effect 
of knock-out mutations or changes in initial conditions in this framework, and the agreement between a real 
system and a Boolean model of it is a strong indication of the robustness of the system to changes in kinetic 
details |Albert and Qthmer, 20031 . 

In discrete models the decision whether a network node (component) will be affected by a synthesis 
or decay process is determined by the state of effector nodes (nodes that interact with it). Typical time- 
dependent Boolean models use synchronous updating rules [Kauffman et al., 2003||Albert and Qthmer, 2003| 
IBodnar, 1997||Kauffman, 1993| , assuming that the time scales of the processes taking place in the system 
are similar. In reality the timescales of transcription, translation, and degradation can vary widely from 
gene to gene and can be anywhere from minutes to hours. Logical models following the formalism intro- 
duced by Rene Thomas [Thomas, 1973[ allow asynchronism by associating two variables to each gene: a 
state variable describing the level of its protein, and an image variable that is the output of the logical rule 
whose inputs are the state variables of effector nodes. Whether the future state variable of a gene equals 
the image or current state variable depends on the update order and, in the absence of temporal information, 
the Thomas formalism focuses on determining the steady states, where the state and image variables co- 
incide IIMendoza et al, 1999| [Sanchez and Thieffry, 2001||Ghysen and Thomas, 2003]|Bemot et al., 2004|. 



The effect of asynchronous updates on the dynamics of the system, however, has not been explored yet. 

In this paper, we present a methodology for testing the robustness of Boolean models with respect to 
stochasticity in the order of updates. Through this, we are also probing the system itself: will individual 
variations lead to unexpected gene expression patterns? In the asynchronous method, the synthesis/decay 
decision is made at different time-points for each node, allowing individual variability in each process' du- 
ration, but more importantly, it allows for decision reversal if the dynamics of effector nodes changes. It 
becomes possible to reproduce, e.g., the overturning of mRNA decay when its transcriptional activator is 
synthesized, a process that synchronous update cannot capture. Thus, replacing synchronous with asyn- 
chronous updates is not merely a technical detail, but rather a fundamental paradigm shift from pointwise in 
time to potentially continuous communication between nodes. Indeed, the effective synthesis or decay time 
for a certain node are determined by the time interval between the latest update of its effector nodes and its 
current update time, and can be any positive fraction of the unit time interval. We propose three algorithms, 
with varying freedom in the relative duration of cellular processes, and find that very short transcription or 
decay times have the potential to derail the wild type development process. 

The steady states of a Boolean model will remain the same regardless of the mechanism of update, but 
its dynamical behavior can be drastically altered due to the stochastic nature of the updates; for instance, 
the same initial state may lead to different steady states or limit cycles. Since the duration of synthesis and 
decay processes is not known, we randomly explore the space of all possible timescales and update orders, 
and derive the probability of different outcomes. Our methods offer a systematic way of exploring generic 
behavior of gene regulatory networks and comparing it to experimentally observed outcomes. To present a 
concrete example, we generalize a previously introduced Boolean model of the Drosophila segment polarity 
genes ||Albert and Qthmer, 2003 j . This model reproduces the wild type steady state pattern of the segment 



2 



WG 





WG 




PTC 


1 


HH 





Figure 1: The network of interactions between the segment polarity genes. The grey background layers il- 
lustrate two neighboring cells, indicating that some interactions in this network are inter-cellular. The shape 
of the nodes indicates whether the con^esponding substances are mRNAs (ellipses) or proteins (rectangles). 
The edges of the network signify either biochemical reactions (e.g. translation,protein interactions) or reg- 
ulatory interactions (e.g. transcriptional activation). The edges are classified as activating (^) or inhibiting 
(H ). Figure adapted from [Albert and Qthmer, 2003| . 

polarity genes as well as the gene patterns of mutants, but its dynamic behavior is not directly comparable 
to that of the real system. Here we show that asynchronous update leads to a much more realistic model that 
gives further insights into the robustness of the gene regulatory network. 



2 The segment polarity gene network in Drosophila 

The Drosophila melanogaster segment polarity genes represent the last step in the hierarchical cascade 
of gene families initiating the segmented body of the fruit fly. While the preceding genes act transiently, 
the segment polarity genes are expressed throughout the life of the fly, and their periodic spatial pattern is 
maintained for at least 3 hours of embryonic development [Wolpert et al., 1998] . The regulatory roles of the 
previously expressed genes such as the pair-rule genes fushi tarazu, runt, even-skipped are incorporated in 
the prepattern (initial state) of the segment polarity genes. The stable maintenance of the segment polarity 
gene expression is due to the interactions between these genes (see FigureQl, and it is a crucial requirement 
in the development and stability of the parasegmental furrows. The best characterized segment polarity 
genes include engrailed (en), wingless (wg), hedgehog (hh), patched (ptc), cubitus interruptus (ci) and 
sloppy paired (sip), encoding for diverse proteins including transcription factors as well as secreted and 
receptor proteins. 

The pair-rule gene sloppy paired (sip) is activated before the segment polarity genes and expressed 
constitutively thereafter [Grossniklaus et al., 1992| Cadigan et al., 1994| . sip encodes two forkhead domain 



transcription factors with similar functions that activate wg transcription and repress en transcription, and 
since they are co-expressed we designate them both SLR The wg gene encodes a glycoprotein that is se- 
creted from the cells that synthesize it |Hooper and Scott, 1992 IPfeiffer & Vincent 19991 . and can bind 



to the Frizzled receptor on neighboring cells, initiating a signaling cascade leading to the transcription of 
engrailed (en) ^Cadigan & Nusse 1997| . EN, the homeodomain-containing product of the en gene, pro- 
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motes the transcription of the hedgehog gene (M) [Tabata et al., 19 92"|. In addition to the homeodomain, 
EN contains a separate repression domain that affects the transcription of ci |Eaton & Romberg 1990| and 
possibly ptc [Hidalgo and Ingham, 1990 Taylor et al., 1993 1. The hedgehog protein (HH) is tethered to the 
cell membrane by a cholesterol linkage that is severed by the dispatched protein, freeing it to bind to the HH 
receptor PTC on a neighboring cell [Ingham & McMahon 2001 J . The intracellular domain of PTC forms a 
complex with smoothened (SMO) in which SMO is inactivated by a post-translational conformation change 
(Ingham 1998). Binding of HH to PTC removes the inhibition of SMO, and activates a pathway that results 
in the modification of CI [Ing ham 1998 |. The CI protein can be converted into one of two transcription 
factors, depending on the PTC-HH interactions. In the absence of HH signaling CI is cleaved to form CIR, 
a transcriptional repressor that represses wg, ptc and hh transcription| Aza-Blanc & Romberg 1999 1. When 
secreted HH binds to PTC and frees SMO, CI is converted to a transcriptional activator, CIA, that promotes 
the transcription of wg and ptc \ Aza-Blanc & Romberg 1999 Ohimeyer & Kalderon 1998] . 

The initial state of the Drosophila segment polarity genes includes two-cell-wide SLP stripes followed 
by two-cell-wide stripes not expressing SLP |Cadigan et al., 1994) , single-cell-wide wg, en and hh stripes 



followed by three cells not expressing them, and three-cell-wide stripes for ci and ptc \ Hooper and Scott, 1992 



[Wolpert et al., 1998J : This pattern is maintained almost unmodified for three hours' (see Fig. |2t), during 
which time the embryo is divided into 14 parasegments by furrows positioned between the the wg and en 
-expressing cells |Hooper and Scott, 1992| . 

The first model of the segment polarity gene network was proposed by von Dassow and collaborators 
Ivon Dassow et al., 2000 1, and is a continuous-state model of 13 equations and 48 unknown kinetic param- 
eters. The main conclusion of the | von Dassow et al., 2000 1 article is that the gene patterns are robust with 
respect to variations in the kinetic constants in the rate laws, thus the essential feature of this network is 
its topology, i.e. the existence and signature (activating or inhibiting) of the interactions. The idea of the 
network topology determining its dynamics was further explored by [Albert and Othmer, 2003| , who used a 
slightly different network reconstruction and assumed synchronous Boolean regulation among nodes. In the 
jAlbert and Othmer, 2003 1 model each mRNA and protein is represented by a node of a network, and the 
state of each node is 1 or 0, according to whether the corresponding substance is present or not. The states of 
the nodes are updated synchronously, and the future state of node i is determined by a Boolean function of 
its current state and the current states of those nodes that have edges incident on it. The updating functions 
are based on the experimental information and on the following dynamical assumptions: (i) the synthesis of 
mRNAs/proteins has the duration of one timestep; (ii) the effect of transcriptional activators and inhibitors 
is never additive, but rather, inhibitors are dominant; (iii) mRNAs decay in one timestep if not transcribed; 
(iv) transcription factors and proteins undergoing post-translational modification decay in one timestep if 
their mRNA is not present; (v) protein-protein binding, such as in the formation of the Patched-Hedgehog 
complex, is assumed to be instantaneous. In summary, the [Albert and Othmer, 2003 1 model assumes that 
gene transcription, protein translation, mRNA and protein decay all happen on a similar timescale, while 
protein complex formation is instantaneous compared to this common timescale. 

The [von Dassow et al., 2000 1 and [ Albert and Othmer, 200 31 models agree in their conclusions regard- 
ing the robustness of the segment polarity gene network. The simplicity of the Boolean rules in the latter 
also allows for the exploration of knock-out mutations and changes in the prepattern of the segment polar- 
ity genes. Starting from the known initial state of en, wg, hh, ptc, ci and SLP, and assuming the null(off) 
state for all other nodes the [Albert and Othmer, 2003 1 model leads to a time-invariant spatial pattern (see 
Fig. 2a) that coincides with the experimentally observed wild-type expression of the segment polarity genes 

'a notable exception includes the refinement of the ptc pattern. 
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Figure 2: a) Top: Illustration of the gene expression pattern of wingless on a gastrulating (stage 9) em- 
bryo. Other segment polarity genes have similar periodic patterns that are maintained for around three 
hours of embryonic development. The parasegmental furrows form at the posterior border of the wg- 
expressing cells |Wolpert et al., 1998| . Bottom: Synthesis of the wild type expression patterns of the seg- 
ment polarity genes (see also text) fHooper and Scott, 1992 Wolpert et al., 1998 |. Left corresponds to an- 
terior and right to posterior in each parasegment. Horizontal rows correspond to the pattern of individual 
nodes - specified at the left side of the row - over two full and two partial parasegments. Each paraseg- 
ment is assumed to be four cells wide. A black (gray) box denotes a node that is (is not) expressed, b) 
Top: wingless expression pattern in an patched knock-out mutant embryo at stage 11 [Tabata et al., 199 2 |. 
The wingless stripes broaden, and secondary furrows appear at the middle of the parasegment, indicat- 
ing a new en-wg boundary. Bottom: Broad striped steady state of the Boolean model, obtained when 
patched is kept off (with the change that ptc and PTC are not expressed), or when wg, en, hh are initi- 
ated in every cell [Albert and Othmer, 2003 1 . This steady state agrees with all experimental observations 
on ptc mutants and heat-shocked genes |Tabata et al, 1992||Gallet et al, 2000| [Martinez- Arias et al, 1988| 
ISchwartz et al, 1995[ jPiNardo et al., T988] [Ingham et al., 1991] IBejsovec and Wieschaus, 1993| |. c) Top: 
wingless expression pattern in an engrailed knock-out mutant embryo at stage 11 [Tabata et al., 19921 . 
The initial periodic pattern is disappearing, and gives rise to a non-segmented, embryonic lethal phe- 
notype. Bottom: Non-segmented steady state of the Boolean model, obtained when wg, en or hh 
are kept off, or cell-to-cell signaling is disrupted [Albert and Othmer, 2003 1 . This steady state agrees 
with all experimental observations on wg, en, hh mutants [Tabata et al., 19921 [DiNardo et al., 1988] 
ISchwartz et al, 1995[ [Hidalgo and Ingham, 1990| [Gallet et al., 2000) . Gene expression images obtained 
from http://www.fruitfly.org (a) and [Tabata et al., 19921 (b,c). 
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during stages 9-11. Indeed, wg and WG are expressed in the most posterior cell of each parasegment, 
while en, EN, hh and HH are expressed in the most anterior cell of each parasegment, as is observed ex- 
perimentally |Ingham 1998 [Tabata et al., 1992] , ptc is expressed in two stripes of cells, one stripe on each 



side of the en-expressing cells, the anterior one coinciding with the wg stripe [Hidalgo and Ingham, 1990 



Hooper and Scott, 1992 1. ci is expressed almost ubiquitously, with the exception of the cells expressing en 



|Eaton & Romberg 1990| . CIA is expressed in the neighbors of the HH-expressing cells, while CIR is ex- 



pressed far from the HH-expressing cells |Aza-Blanc & Romberg 1999 1. The model indicates that knock- 



out mutations in en, wg, hh cause the non-segmented gene pattem shown on Fig. 2b, which agrees with 
experimental observations. Indeed, the hh expression in en null embryos starts normally, but disappears 
before stage 10 ^Tabata et al., 1992J . In wg null embryos, en is initiated normally but fades away by stage 
9, as observed by DiNardo et al. (1988), while ci is ubiquitously expressed [Schwartz et al., 1995| . In hh 



mutant embryos the wg expression disappears by stage 10 [Hidalgo and Ingham, 1990 1, as does the expres- 



sion of ptc, and there is no segmentation [Gallet et al., 2000| . On the other hand, ptc knockout mutations 
or overexpressed initial states lead to the broad-striped pattem of Fig. 2c^. Indeed, experimental results 
indicate broad en, wg and hh stripes [ Tabata et al., 19921 [Gallet et al., 2000| [Martinez-Arias et al, 19881 
and Gallet et al. (2000) find that a new ectopic groove forms at the second en — wg interface at the 
middle of the parasegment. Also, ci is not expressed at this ectopic groove ISchwartz et al., 1995| . In 
heat-shock experiments the wg and ptc stripes expand anteriorly when hh or en are ubiquitously induced 
IGallet et al., 2000| , while narrower ci stripes emerge after a transient decay of ci [Schwartz et al., 1995| . 
Intriguingly, the [Albe rt and Qthmer, 2003 1 [ model finds that a knock-out mutation of ci does not change the 
en, wg, hh patterns but disrupts ptc expression; experiments indicate that the segmental grooves are present 
and wg is expressed until stage 11, but ptc expression decays [.Gallet et al., 2000 [. In summary, the simple 
synchronous Boolean model [Albert and Qthmer, 2003 [ captures perfectly the wild type and mutant expres- 
sion patterns of the segment polarity genes, and thus serves as a good starting point for a more realistic 
model that relaxes the assumption of synchronicity. 

We focus our attention on a single parasegment of four cells, thus the total number of nodes we con- 
sider is 4 X 13 = 52. We use the same interaction topology and logical rules as the synchronous model 
[Albert and Qthmer, 2003| , but instead of assuming that the states of all nodes are updated simultaneously, 
we update the state of each node individually (see TableQl- To maintain the highest generality, we incorpo- 
rate possible cell to cell variations in synthesis and decay processes.^ 

Throught the text, the notation ''wg\" or "wgi{t)" represent the state of wingless mRNA in the first 
cell of the parasegment at time t. Similar notations apply for other mRNAs and proteins. There are 4 cells 
in each parasegment, and we adopted periodic boundary conditions, meaning that: no(ie4+i = nodei and 
nodei-i = node^. The wild type initial state corresponds to: 

wgl = l, en? = l, hh\ = l, ptc^,3,4 = 1, ci^2,ZA = '^ (1) 

and the remaining nodes are zero. The asynchronous model represented in Tabled exhibits the same steady 
states as the synchronous model developed in [Albert and Qthmer, 2003| . Note that three of the four main 
steady states agree perfectly with experimentally observed states corresponding to wild type, en, wg or hh 
mutant and ptc mutant embryonic pattems [Tabata et al., 19921 [DiNardo et al., 1988] [Schwartz et al., 1995| 

^The only difference between the ptc mutant and heat-shock pattern is that the former does not express ptc and PTC 

^We follow the [Albert and Qthmer, 2003 J model in assuming very short timescales for PTC-HH binding and SMO activation, 

and consequently in Figure Q and in the regulatory rules we connect the CI posttranslational modifications to HH signaling. We 

have verified that this assumption can be relaxed without any qualitative changes in the results. 
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Table 1 : Regulatory functions governing the states of segment polarity gene products in the model . Each 
node is labeled by its biochemical symbol and subscripts signify cell number. The times signify the last 
time node j was updated before t. 

Node Boolean updating function in the asynchronous algorithm 



SLPi SLPi{t) 



if iG{l,2} 

1 if i G {3, 4} 

WQi wgi{t) = {CIAiijciA) and SLPi{TsLp) and not CIRiijciR)) 

or [wgi{Ty,g) and {CIAi{TciA) or SLPi{TsLp)) and not CIRiiTcm)] 
WGi WGi{t) = wQiiT^g) 

erii emit) = {WGi-i{TwGi) or WGi+i{TwG2)) and not SLPi{TsLp) 
ENi ENiit) 

hhi hhi{t) = ENi{TEN) and not CIRi{TciR) 
HHi HHi{t) = hhi{Thh) 

ptci ptci{t) = GIAi{TciA) and not ENi{TEN) and not CIRi{TciR) 

PTGi PTGiit) = ptCiiTptc) or {PTGi{TpTc) and not HHi_i{THm) and not HHi+i{THH2)) 

cii cii{t) = not ENi^TEN) 

Gli Chit) = cii{Tci) 

GIAi CIAiit) = Cliirci) and [not PTGi{TpTc) or HHi^i{THHi) 

or HHi+i{THH2) or hhi^i{Thhi) or hhi+i{Thh2)] 
GIRi CIRiit) = Gh{Tci) and PTG^{TpTc) and not HHi_i{THm) and not HHi+i{THH2) 

and not hhi-i{Thhi) and not hhi+i{Thh2) 



Hidalgo and Ingham, 1990l|Gallet et al, 2000| [Martinez- Arias et al, 1988||Bejsovec and Wieschaus, 1993 
Ingham et al., 199l[ Hooper and Scott, 1992 Wolpert et al., 19981 . ^ summary is presented in Tabled 



3 Randomly perturbed timescales 

As in the context of parallel computation systems, the fundamental difference between synchronous and 
asynchronous updates is at the level of task coordination and data communication among nodes in a net- 
work [Bertsekas and Tsitsiklis, 1989] . Synchronous algorithms are highly coordinated: at pre-determined 
instants, all the nodes "stop" and exchange the current information among themselves. For instance, sup- 
pose there are N nodes, where each node i "computes" the state of variable x-i, according to a function 
fi{xi,X2, . . . ,xn) (i = 1, • • • , N). When all the nodes have finished phase k, they exchange their 
current states, x^, and then proceed to phase k + 1, that is 

fc+l _ f.(^k k k \ 

Asynchronous algorithms, on the other hand, admit a greater flexibility at the level of process coordination. 
Each node is allowed to have its own "computation rate", that is, during any time interval [ta^ifo]. node 
i may be updated only once, while node j may be updated £ > 1 times. In this case, communication 
delays between nodes may occur, and some possibly outdated information may be used: for instance, node 
j uses the same value Xi{ta) throughout its i updates in the interval [tai^fe]- However, an overall gain in 
efficiency in achieving the final result may be expected. For instance, in our example, the wild type steady 
state is reached in less than 4 steps with the asynchronous algorithms (see Sections |4j [Sj, while with the 
synchronous algorithm 6 steps are needed ^Albert and Othmer, 2003| . 
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Table 2: Complete characterization of the model's steady states. 



Steady state 


Expressed nodes 


wild type 


wg4^, WG/j^, eni, ENi, hhi, HHi, 

ptC2,4, PTC2,3,4' Ci2,3,4, Cl2,3,4, CIA2,A, CIR3 


broad stripes 


ptC3^4, PrC3,4, d3,4, C/3,4, CIA^^i 


no segmentation 


cn,2,3,4, C'/i,2,3,4, -PrCi, 2,3,4, C'/i?i,2,3,4 


wild type variant 


WQi, WG4, eni, ENi, hhi, HHi, 

ptC2,4, -PrCi,2,3,4, Ci2,3,4, C'/2,3,4, ^1^2,4, CIR3 


ectopic 


wg3, WG3, en2, EN2, /1/12, HH2, 

ptCi,3, PTCi,3A' Cii,3,4, C/1,3,4, C/Ai,3, CIRa 


ectopic variant 


wg3, WG3, en2, EN2, /1/12, HH2, 

ptCi,3, -P7'C'i,2,3,4, Cii,3,4, C/1,3,4, CMi,3, CIRi 



In general, we may say that node i updates its state at times: 

Tl,T^,...,Tt,... fee No, 
and the local variables, Xi, are updated according to: 

Xi[Tl^] = n{x,[T^,],...,XN[r'^,]), (2) 

where t^, is defined as 

rj^j = the latest available communication to node i, from node j. 

There is usually a distinction between totally or partially asynchronous algorithms: the latter impose an 
updating constraint (every variable is updated at least once in any interval of a fixed length), while the 
former simply ensure that a variable is updated infinitely many times. 

In a first numerical experiment we consider a totally asynchronous algorithm, with the highest degree 
of individual variability in each process' duration. The time unit of the synchronous model is randomly 
perturbed, so that the set of updating times for each node i (1 < i < A'^) is of the form 

Xk+^ =T^^ + l + erf, ken, 

where are random numbers generated at each iteration, out of a uniform distribution in the interval [—1,1]. 
The value e S [0, 1) is the magnitude of the perturbation (the case e = coincides with the synchronous 
algorithm). At any given time t, the next node(s) to be updated is(are) j such that Tj = mmi^k{T^ > t}, 
for some £. Since the duration of synthesis and decay processes is not known, through this algorithm one 
may randomly explore the space of all possible timescales and update orders, and derive the probability of 
different outcomes. The set of updating times {Tf^, /c S No} may vary with each execution of the algorithm, 
so an element of stochasticity is naturally introduced. 

Always starting from the wild type initial condition Q, this experiment was conducted over a wide 
range of perturbations (10~^^ < e < 0.65), and 30000 trials were executed for each e. The results (see Fig- 
ure ISj show that all of the model's steady states may occur with a certain frequency: the wild type pattern 
with only 57%, followed by the broad-striped pattern (24%) observed in heat-shock experiments and ptc 
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Figure 3: Fragility of the regulatory network. With the totally asynchronous algorithm the wild type initial 
state can lead to one of the six distinct steady states. Each * corresponds to an e perturbation of the unit 
time-step. Note that the e limit does not give the same results as a synchronous update, demonstrating 
the fundamental difference between synchronous and asynchronous models. 

mutants fGallet et al., 2000] and by the pattern with no segmentation (15%) observed in en, hh or wg mu- 
tants [Tabata et al, 1992J , the latter two corresponding to embryonic lethal phenotypes IGallet et al., 2000) . 



We observe that each of the steady state patterns occurs with a frequency which is independent of the 
value of e, for e < 0.15. This may indicate that it is the order in which the protein and mRNA nodes 
are updated that determines the steady state pattern. In order to test this hypothesis, we designed a second 
experiment assuming that 

(Al) Every node is updated exactly once during each unit time interval {k, + 1] (A; = 0, 1, 2, . . .), accord- 
ing to a given order (j)''. 

This order cp^ is a permutation of {1, . . . N}, chosen randomly (again out of a uniform distribution over the 
set of all A^! possible permutations) at the beginning of the time unit k. Then we have 

T^^ = N{k- I) + (i), ken, 

so that < implies Tj^ < Tf^, and node j is updated before node i. The partially asynchronous 

algorithm leads to the same patterns, with incidence rates very similar to those observed with the totally 
asynchronous algorithm (see Table |3ll. 

These results indicate the fragility of the wild type gene pattern with respect to changes in the timescales 
of synthesis and decay processes. While more than half of the random timescale combinations still lead to 
the expected outcome, a considerable percentage results in loss of the prepattern and an inviable final state. 



3.1 Imbalance between CIA and CIR 

Further analysis shows that the divergence from wild type can be attributed to an imbalance between the 
two opposing Cubitus Interruptus transcription factors (CIA, CIR) in the posterior half of the parasegment. 
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Table 3: The frequencies of the six steady states observed in the partially asynchronous model confirm those 
observed for the totally asynchronous model. The frequencies are computed from 30000 executions. 



Steady State 



Incidence 



wild type 
broad stripes 
no segmentation 
wild type variant 
ectopic 

ectopic variant 



56% 

24% 

15% 

4.2% 

0.98% 

0.68% 



Indeed, the expression of CIA and CIR in both the broad stripes and the no segmentation patterns is clearly 
distinct from that in the wild type pattern. In the next set of numerical experiments, we explore the effects 
of CIA/CIR expression in the formation of the final pattern. 

In wild type, the two Cubitus Interruptus proteins, CIA and CIR, are expressed in different cells of the 
posterior part of the parasegments, namely. 



and the maintenance of these complementary ON/OFF states is essential in the wild type pattern. To investi- 
gate the effect of an imbalance between the two Cubitus Interruptus proteins, we considered two disruptive 
cases: the (transient) overexpression of CIR, or the (transient) overexpression of CIA and absence of CIR in 
both posterior cells. 

More precisely, in the totally asynchronous algorithm (choosing e = 0.1), we transiently imposed an 
expression pattern for the Cubitus proteins as follows: 



where r is the duration of the transient. The overexpression starts after three unit time steps. The duration 
of the transient was: 



so when r = the results of the general totally asynchronous algorithm are recovered. 

Our results show that even a small transient imbalance between CIA and CIR causes a clear bias towards 
a mutant state: the broad stripes mutant in case (a), or the no segmentation mutant in case (b). Thus any 
perturbation that leads to such an imbalance has as severe effects as a mutation in ptc (causing the broad 
striped pattern) or either of en, wg or hh (causing the nonsegmented pattern). 

These numerical experiments also open the way to many other questions: are there particular sequences 
that lead to a given steady state? How is the evolution from the initial to steady state? How robust is the 
asynchronous model with respect to initial conditions? 



CM3 = 0, CIAi = 1, 
C/i?3 = 1, CIRi = 0, 



(a) C/^l 4 = 1 and C/i?| 4 = 0, for t G [3, 3 + r]; 

(b) C/i?|_4 = 1, for i € [3, 3 + r]. 



r G {0,0.3,0.75,1.5,2.75,3}, 
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Figure 4: Bias towards mutant states. The x-axis represents the duration of the transient, r (in unit time 
steps). The incidence probabihties were computed over 20000 trials, (a) The case CIA\ 4 = 1 and 
CIR\ 4 = leads to the broad striped pattern, (b) The case CIR\ 4 = 1 leads to the no segmentation 
pattern. 

4 Timescale separation uncovers robustness of the model 

In both of the previous algorithms we assumed no bias towards a preferred protein/mRNA updating sequence 
and, as a result, an unrealistic divergence from the wild type pattern is observed, with high incidence of 
inviable states. Based on the fact that post-translational processes such as protein conformational changes 
or complex formation usually have shorter durations than transcription, translation or mRNA decay, we 
introduce a distinct timescale separation by choosing to update proteins first and mRNAs later. This leads to 
a model which is very robust, in the sense that the wild type pattern occurs with a frequency of 87.5% and 
only one other steady state is observed, the broad striped pattern, with a frequency of 12.5%. We completely 
characterize this model by theoretically showing that only two of the six steady states are possible (and occur 
with well determined frequencies), and identifying the order of updates that leads to divergence from wild 
type. We also show that the wild type state is really an attractor for the system, while the pathway to the 
broad stripes state may show oscillatory cycles. 
Assuming that 

(A2) All the proteins are updated before all the genes, 

the fc-th iteration of the two-timescale algorithm proceeds as follows: 

(A3) At the begining of the /c-th time unit, generate a random permutation, of {1, . . . L}, and a random 
permutation, 0,^^^^ of {L + 1, . . . N} (using a uniform d istribution over, respectively, the sets of 
L\ and {N — L + 1)! possible permutations). Then the nodes are updated in the order given by 

<l^^ = {(Ppron (ttRNA)' according to ©, with 




As an example, suppose that 

N = 5, L = 3, 0L = {2,1,3}, <pl,, = {5,4}- 
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Table 4: Regulatory functions governing the states of segment polarity gene products in the two-timescale 
asynchronous algorithm. Each node is labeled by its biochemical symbol, subscripts signify cell number and 
superscripts signify timestep. Although the updating time of each node varies, each function can be written 
by using the states of effector nodes at the previous or current timesteps. The individual times ti . . . tiQ can 
take the values {t — l,t}. 

Node Boolean updating function in the two-timescale algorithm 



wgl = {CIA\ and SLPl and not CIRf) or [wg^r^ and {CIA\ or SLPl) and not C/i?*] 

WGi WG\ = wg^-^ 

em en\ = {WG\_^ or WG\^^) and not SLPl 

ENi ENl = en*-^ 

hhi hh\ = ENj and not GIR\ 

HHi HHl = hhl-^ 

ptci ptc\ = GIA\ and not ENf and not GIR\ 

PTCi PTGj = ptcl'^ or {PTGj-^ and not HHI'_^ and not HH^l^) 

oil ci\ = not ENl 

Gh Gil = ci\-^ 

GIAi GIA\ = Gll^ and (not PTC*^ or HHI'L^ or HHII^ or /i/i*^}) 

GIR^ GIR*, = Gil' and PTC** and not HHI'' , and not and not /i/i*7 j 

Then, (jr^ = {2, 1, 3, 5, 4}, and = 2, = 1, Tl = 3, T} = 5, = 4. The nodes are updated as follows 
(for simplicity of notation, we will write := XilT^]): 



Some general inferences about the updating rules can be made. For example, the translation process only 
depends on the presence of the transcript, which is decided in the previous time unit, thus Prot^ = mRNA^^^. 
The beginning of a transcription process depends on the presence of transcription factors, and since mRNAs 
are updated after proteins, mRNA* = Prot*. The outcome of post-translational processes depends on the 
order of updates, for example the rule for a binding process will be Complex* = Prot*^ and Proi^, where t\ 
and t2 can be either t — 1 or t (see Table HJ. 

4.1 Two steady states 

The trajectory of the system is thus defined by a sequence of permutations (obtained as described in A3) and 
the corresponding sequence of states: 



We will show that for pre-pattems that satisfy wg\ = 1, cvl = and ptcl = (which include the pattern 
observed in the wild type at stage 8), the only possible steady states for system Q are the wild type pattern 
experimentally observed at stages 9-11, and a mutant with broad wg stripes. We assume that all the proteins 
are absent initially (at T = 0), and that the which sloppy pair gene is maintained at a constant value: 




fl{Xi, X2, X3, X^, X5), 
fs{x\,X2,X^,X^, X5), 
f5{x\,X2,x\,X^, X5), 
f4{x\,X2,x\,X^, X5). 



{((>''}, {x'^} for A: = 0,1, 2,.... 



(3) 
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SLPi 2 = and SLP^^^ = 1 This pattern for SLP is responsible for permanent absence (or expression) 
of some of the segment polarity genes, and corresponding proteins, in certain cells of the parasegment. By 
direct inspection of the model, it follows that 

wgl2 = 0, wgl2 = 0, WGl2 = 0, for T > (4) 

en|^4 = 0, ENJ^^ = 0, for T > 0, (5) 
/i/i|^4 = 0, HHJ^^ = 0, for T > (6) 



and 



^ = 1, ^ d|^4 = 1, for T > 0, and Cll^ = l,for T > 1 (7) 
0, ^ ciL = 1, for T > 1, and CiL = l,for T > 2. (8) 



ci 



3 4 — u, ^ ct3_4 — ±, lui -t ^ ±, aiiu <-v-i34 — ±, 

The next statement reflects the fact that the effect of wg^ activating eni propagates to inhibit cii which then 
eliminates all forms of CI from the first cell. 

Fact 1. Assume that wgl = 1, ci^ = and ptc^ = 0. For any T > 0, if wgl = 1 for all < t < T, then 
C/^ = for all 3 < t < r + 3 and CIR\ = for all < t < T + 3. 

Proposition 4.1. Assume wg^ = 1, ci^ = and ptci = 0. Under assumptions A1-A2, wgj = 1, for all 
T > 0. 

Proof. We will argue by contradiction. Suppose that there do exist times t > 1 with wg\ = 0, and let T be 
the minimum of such times, that is, 

wgJ = and wg^ = 1, for all < i < T. 

^From the model's equations, together with assumptions A1-A2: 

WGi = wg\-\ 
wgl = {CIA\ and not CIR\) or {wg*^^ and not CIR\), 

for all t > 1. So, it follows that 

WG\ = 1, forallO<t<r, (9) 
CIR\ = 0, for ah < i < T, and CIRJ = 1. (10) 

Now, from Fact^it also follows that 

CIR{ = for aU < t < T + 2. (11) 

The equation for CIR4 is: 

CIRI = Cll" and not [not PTCl'' or HHI'' or HHI" or hhl'^ or /i/i*"^] (12) 

(13) 

(where ta, ■ ■ ■ ,td £ {t,t — 1} depend on the permutation (/>*). Recall also that 

hh\ = E Nl and note I R\ (14) 
ENl = en\-^ (15) 
or VFG* 



en\-^ = WGi-^ or WGi-^. (16) 
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From (fT2l : 



CIR^ = 1 ^ hhj-^ = 0, 

and then from dTTT l and (fTlT) : 

hhj'^ = ^ EiVj^-^ = 0. 

Now by equations (IT5lfT6b : 

EN^-^ = 0, ^ enf "2 = ^ VFGJ"^ ^ g and T^G^'^ = 0, 
which contradicts equation Thus, it must be that wgj = 1 for all times T, as we wanted to show. □ 
The following are now immediate conclusions from the model. 

Corollary 4.2. CIRJ = for all T > 0, enf = 1 and WGj = 1 for all T > 1. EN{ = 1, cif = and 
hhj = 1 for all T > 2. Clf = and HHf = 1 for all T > 3. And finally, CIAj = CIRJ = for all 
T > 4. 

Corollary 4.3. ptcf = and PTCf = for all T > 0, and Cli?^ = for all T > 3. 

In conclusion, from Proposition 14.11 it is clear that neither the no segmentation nor the two ectopic 
patterns are steady states of the system Q under assumptions A1-A2, because all of these states imply 
wg4 = 0. In addition. Corollary 14. 31 shows that the wild type variant, where PTC is ubiquitous, cannot be 
a steady state. Also, any of the states with wgi 2 = 1 is immediately prevented by the initial condition (|4]i. 
This leaves only the "regular" wild type or the mutant with broad wg stripes. 



4.2 Divergence from wild type 

Under assumptions A1-A2, divergence from the wild type pattern occurs if and only if the first permutation 
(in particular (pj^^^) is of a particular form. Thus, convergence (or divergence) to the wild type pattern is 
decided at the first iterate (T = 1). 

Recall that the wild type pattern requires wingless not to be expressed in the third cell (wg^ = 0). The 
next Fact (proved in the Appendix) essentially says that a stable wg^ = induces the absence of both 
engrailed, hedgehog in the second cell, as well as the absence of CIA3, and maintains the expression of 
PTC3. 

Fact 2. Assume pic!] = 1 and = 0. 
(a) Let r > 1. If wgl = for all < t < T, then 



0, EN^ = 0, 0<t<T + 2, 

/i/i2 = 0, 1 < i < T + 2, 

HH*2 = 0' 1 < t < T + 3, 

PTCl = 1, 1 < t < r + 3, and 

CIAl = 0, CIRi = 1, 2 < t < r + 3. 

(b) Furthermore, if ci^ = 0, then also CIAl = and part (a) holds for any T > 0. 
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With the help of this Fact, we establish that wg^ may become expressed only at the first iterate or else 
it is never expressed. Thus the two timescale model provides a strong natural restriction on the formation 
of an inviable state: if wg^ = 0, then wg'^ = for all T > 0, implying that such trajectories will never 
converge to the broad striped pattern. 

Proposition 4.4. Assume that the initial condition satisfies wg^ = 0, ptcl = 1, hhl 4 = 0, and d!] = 1. 
Then wgj^ = 1 and wg'^ = for all < T < Ti, only if Ti = 1. 

Proof. To obtain wgj = 1 with wgj^^ = it is necessary that 

wgj = CIAJ and not CIRJ =^ CIAJ = 1 and CIRI = 0. 

But, if wg'^ = for T = 1, then, by FactEJ the activator CIA^ is zero for T = 2, 3, 4. Then (by induction 
on T) expression of wg'i is prevented at any later time. □ 

In addition, it is possible to completely characterize the updating permutation {(j)^) that leads to wg\ = 1 
and, as a consequence, exactly compute the probability of divergence (hence convergence) to the wild type 
steady state (Section l43b . 

Proposition 4.5. Assume that assumptionts Al and A2 hold. Assume that the initial condition satisfies 

wgl = 0, ptcl = 1 and hh\^ = 0. 

(a) If cil = 0, then wgl = 0. 

(b) If d!| = l,then wg\ = lif and only if the permutation cf)^ satisfies the following sequence among the 
proteins CI, CIA, CIR and PTC: 

CIR3 Ch CIA3 PTCs, 

Ch CIR3 CI A3 PTC3, (17) 

CI3 CI A3 CIR3 PTC3, 

while the other proteins may appear in any of the remaining slots. 

Proof. Part (a) follows immediately from Fact|2lb). To prove part (b), we start by noticing that, because 

SLP3 = 1 and wg^ = 0, 

wgl = CIAl and not CIRI, 

so that 

wgl = l ^ CIAl = 1 and CIRI = 0. 

Following assumptions A1-A2, the model's equations for CIA^ and CIR\ are given by: 

CIA\ = Cll^ and [not PTC*" or HH\- or HHI^ or hh^ or hh^] , 
CIRI = CIl'' and not [not PTCg" or HH'^'' or HHl" or hh^ or hhl] , 

where ta, ■ ■ ■ ,Sa G {0, 1} and depend on the permutation 0^. These expressions may be simplified by 
observing that: (a) 4 = 0, and thus also (b) HHI^'I = 0. Therefore, 

CIAl = Cll^ and not PTC^", 
CIRl = CI^- wdPTCl\ 

The values for Cf/ and PTCg'^ are determined by: 
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1. CI^ = and Cli, = cf^ = 1, 

2. PTC^ = and PTCl = ptc^ or [•••] = 1, since ptc^ = 1, 

and recall that both CIA^ = and CIR^ = 0. Therefore, it is necessary that CIs is updated before CIA^ 
and PTCs is updated after CIA^ , because otherwise CIA^ = 0. Finally, CIR^ must be updated before 
PTC3, because otherwise CIRl = 1. In other words: 

ta = 1, h = 0, Sa G {0,1}, Sb = 0. 

It is easy to see that any of the sequences ( I17t is also sufficient to obtain wgl = 1. □ 

Finally, we will show that whenever wgs becomes expressed at time T = 1, it is afterwards periodically 
expressed, every third step. Such trajectories cannot converge to the wild type pattern. In other words, initial 
permutations of the form (fTTl are not included in the basin of attraction of the wild type pattern. 

Proposition 4.6. Assume wg^ = 1, cil = and ptcl = 0. For any T > 1, if wg^ = 1, then wg^~^^ = 1. 

Proof Recall that, by Proposition gl] wg^ = 1 for all t > 0. By Corollary O CIi?| = for all t > 3. 
Now, pick any T > 1 and assume that wgj = 1. Then 

where the last implication follows from Fact|3] Then (using either Q or ©) 

CIAl+'-^ = not PTC3" or HH^^ or hhl+^ , 
CIRJ^^ = PTCI; and not HH'^" and not 

(where ta, tb, Sa, s;, G {T + 2, T + 3}, and depend on the permutation (f)'^^^). So hh^'^'^ = 1 implies 

CIAl+^ = 1 and CIRI+^ = 0, 

and therefore wg'^'^^ = 1, as we wanted to show. □ 

Whenever wg^ is not expressed, some other nodes also stabilize, after the appropriate number of itera- 
tions. These are summarized next. 

Corollary 4.7. Assume that wg^ = for all t > 0. Then WG^ = 0, en^ = for all t > 1. EN^ = 0, 
hhl = and ci2 = 1 for all t > 2. Finally, HH^^ = and Ch = 1 for all t > 3. 



4.3 Probability of convergence to wild type 

The wild type pattern is in fact an attractor for the asynchronous model: every trajectory which is not of the 
form (fnt converges to the wild type pattern (see AppendixlBjl. The probability that this happens is therefore 
determined by counting all the possible states of the form ( fTTt : 

# permutations as in ( fTTt 

Prob(wild type) =1 

Total # permutations 

Let L be the number of protein nodes to be updated at each iterate (there are 9 proteins in each of the four 
cells so L = 36). Out of the L proteins, only 4 need to satisfy one of particular sequences (fTTb in their 
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relative positions. So let be the number of possible permutations satisfying any of the sequences (flTl . 
Then 

Ml (L-4)! 



Prob(wild type) = 1 



L! 



The next Proposition is proved in the Appendix and shows that, in fact, this is a constant number, indepen- 
dent of L. That is. 



Proposition 4.8. For any L > 4, 

3! 



Prob(wild type) = 0.875. 



2 ^ \ 3 / 2 . ^ 

P=4 \ / J=l 



and 



(L-4)! ^ 1 
L\ 8' 



5 A Markov chain process 

As a Boolean model, there are only a finite set, say S, of distinct states (in the total state space {0, 1}^) 
reachable by the system. Starting from any state Sa € S, each permutation of {1, ... , N} takes the system 
to some other state 5{, € S. It is possible to theoretically identify all the distinct intermediate and final states 
of the system as well as all the possible transitions after one iteration. Thus the asynchronous algorithm 
consisting of the N node functions Q together with assumptions (Al), (A2) and (A3) may be characterized 
as a Markov chain process, by identifying the d distinct states 

S = {Si,S2, ■ ■ ■ ,Sd}, 

and the d x d transition matrix P, where 

Pij = probability of a transition from state Si to state Sj . 

The probabilities Pij are simply the fraction of the total number of permutations that (in one iterate) trans- 
form the state Si into state Sj. The matrix P is a stochastic matrix, since all its rows add up to 1. A state 
Sa with the property that all permutations leave the state unchanged (that is, Paa = 1 and Paj = for 
j 7^ 0) is called an absorption state of the Markov chain, and it is also a steady-state of system In the 
asynchronous model there are only two absorption states, corresponding to the wild type and broad wingless 
stripe mutant patterns, as described above. Isolating the two rows and columns that correspond to these 
absorption states, the transition matrix may be partitioned as 



P 



Pa 
Ra P 



where P is of size {d — 2) x {d — 2). It is a well know result(see any standard book on probability theory, 
for instance [Feller, 1970| ) that / — P is an invertible matrix, and 



Ti 

T2 



(J-P) 



-1 



d-2 
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or Tj = 1 + J2'j=i Pij Tj. The values Tj provide an estimated time for absorption when the chain starts 
from state Si (if a is an absorption state, then Ta = 0). In Figure |5l a schematic diagram of transitions is 
shown, together with probabilities and estimated times for absorption. This diagram was obtained from a 
simulation starting with the initial wild type pattern (observed at stage 8 of the embrionic development), and 
following the assumptions (A1)-(A3), as well as the additional 

(A4) Each protein or gene is updated simultaneously in the four cells, 

meaning that there is no cell-to-cell variation in the duration of molecular processes. In this case, the total 
number of possible permutations is a manageable 7! x 5! = 604800: 

$ = {0 = (0^^^,, </,„,^_,) : 0^,.^, is a permutation of {1, 2, 3, 4, 5, 6, 7} , 

(l)„RNA is a permutation of {1, 2, 3, 4, 5} } 

The total number of distinct states of the Markov chain (under assumptions (A1)-(A4)) is d = 48. The 
transition probabilities matrix was computed exactly, by counting all the 7! x 5! transitions from each of the 
48 states. 

It is clear from this diagram that the decision between convergence to the wild type or the mutant 
patterns is indeed decided at the first iteration, in agreement with Propositions 14.11 and 14.41 Furthermore, 
the diagram shows the possibility of periodic oscillations (of period at least three) in the mutant branch 
(see also Proposition 14.61) . Although in practice the probability of a limit cycle is very small, this prevents 
(theoretical) convergence to the mutant state, and considerably increases the absorption times to the mutant 
state. The robustness of the two timescale model is illustrated by the fast convergence to the wild type pattern 
(expected time to convergence is 4 steps), contrasted with the long and oscillation-strewn path toward the 
broad striped pattern (expected time to convergence is 15 steps). 



6 Identifying minimal pre-patterns 

A necessary condition for convergence to the wild type is thatptc3 = 1. Otherwise the trajectory immedi- 
ately fails to enter the basin of attraction of the will type state: 

Fact 3. Assume ptc^ = 0. Then wgj = 1, for some T € {1, 2, 3}. 

Proof. Note that ptc'^ = implies PTCg = 0. Using Q and Q, ©, the equations for ptc^ and CI A3 
simplify to 

CIAl = not CIRl = not PTCg" or ifi?*' or hh'f\ t > r 

ptcl = CIA\ and not CIR\ = CIA\, t > 1 
PTCl = ptcl'^ or [PTC^f^ and not iJF*" and not hh^-^, t>l, 

where r = 2 (respectively, r = 3), if ci'^ = 1 (respectively, ciSj = 0). Consider first the case ci^ = 1. The 
activator protein will be turned on either at the first or second iterations: CI A3 may become activated at 
t = 1 because PTC^ = 0, (if the permutation (p^ is such that CI3 is updated before CIA^). If CIA^ is not 
actived at t = 1, then it certainly is actived at t = 2 (because CIAj^ = implies ptc^ = and PTC^ = 0). 

1 2 

A similar argument shows that CIR^ = 0. 

Consider next the case 0^3 = 0: we have CIA\ = CIR\ = 0, and CIA^ is turned on at the second or 
third iterations, by a similar argument as before. Therefore, the wingless gene is also expressed after CIA^, 
at r = 1,2 or 3. □ 
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MUTANT 



Figure 5: Robustness of the regulatory network modeled with the two-timescale algorithm. There are 48 
states reachable from the wild type initial state. The arrows are labeled by the transition probabiUties be- 
tween states (if unlabeled, the probabiUties are 1), and the expected times to absorption into the correspond- 
ing steady state are indicated between square brackets. 
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Another necessary condition for convergence to wild type, is that 

wg^^ = 1 or en\ = 1 or ci\ = 1. 

Otherwise, the trajectories cannot converge to the wild type nor to the mutant steady states. In this case, the 
only possible steady state is a "lethal" state, where expression of PTC, ci, CIA and CIR is ubiqitous and 
all others are absent. 

Proposition 6.1. Assume wg^^ = 0, en^ = 0, ptc^ = 1 and ci^ 4 = 0. Then 

(a) wg-i = 0, for all t > 0. 

(b) PTCi = 1 for all t > 4. 

(c) wg4 = for infinitely many t. 

Proof. Part (a) follows from Fact|3b), for any T > 0: 

wgl = 0, t<T CI A3 = 0, < t < T + 3. 

Thus, by induction, wgl = for t > 0. 

Part (a) and Factglimply hh^ = and HH^2 = for all i > 0, so that 

ptc{ = CIA\ and not CIR{ and not EN^ 
PTC{ = ptc{-^ or PTCl-^, t>l 
CIA\ = C/{" and not P^C{^ t>l. 

If ptc? = 1, then it is clear that PTCl = 1 for all t > 1. Consider the case ptc^ = 0. Then PTCi = as 
long as ptci = 0. Note that ci^ = implies C/j = 0, CIA\ = and also wg\ = 0. Then wg^^"^ = and 
en°'^ = imply STV^ '^'^ = and CI^'^ = 1. So, either at T = 2 or T = 3, we will have CIAj = 1, and 
therefore also ptc[ = 1 and PTCf^^ = 1. Thus, PTCl = 1 for t > 4, proving part (b). 

Finally, to argure by contradiction, suppose that wg\ = 1 for t > Tq. Then by CoroUarv 14.31 PTCl — ^ 
for t >Ti) > Ta, which contradicts part (b). Hence, wg^ cannot become permanently on. □ 

By Proposition 14. II together with Fact|3b), a sufficient condition for convergence to wild type is 

wgl = l, pte|] = l, pte? = 0, 0^3 = 0. 
Another sufficient condition (which allows the presence of cubitus in the third cell) is 

wgl = 1, ptcl = 1, PTCl = 1. 

The argument in the proof of Proposition 14.51 shows that, if PTC^ = 1 then wg^ = 0. Then, by Proposi- 
tion l4.41 it follows that wg^ = for all times. 

In conclusion, while the wild type initial state allows for an ambiguity in the final states, we find that a 
remarkably minimal prepattern, consisting of 'wg4 and ptc^, is sufficient to guarantee the convergence to the 
wild type steady state. In other words, the initiation of two genes in two cells is enough to compensate for 
initiation delays in any and all other genes, irrespectively of the variations in individual synthesis and decay 
processes. This suggests a remarkable error correcting ability of the segment polarity gene control network. 
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7 Conclusions 



In summary, we proposed an intuitive and practical way of introducing stochasticity in qualitative mod- 
els of gene regulation. We explored three possible ways of incorporating the variability of transcription, 
translation, post-translational modification and decay processes (see Table |5] for a comparison between the 
synchronous and three asynchronous algorithms). Applying our methods on a previously introduced model 
of the Drosophila segment polarity genes gave us new insights into the dynamics and function of the interac- 
tions among the segment polarity genes and, through it, into the robustness of the embryonic segmentation 
process. Our results suggest that unrestricted variability in synthesis/decay/transformation timescales can 
lead to a divergence from the wild type development process, with an expected divergence probability of 
45%. On the other hand, if the duration of post-translational transformations is consistently less than the du- 
ration of transcription, translation and mRNA/protein half-lives, the wild type steady state will be achieved 
with a high probability, despite significant variability in individual process durations. We find that a remark- 
ably sparse prepattem is sufficient to ensure the convergence to the wild type steady state of these genes. 
This dual behavior, robustness to changes in the initial state but fragility with respect to temporal vari- 
ability, is reminiscent of Highly Optimized Tolerance, a feature of highly structured, non-generic complex 
systems with robust, yet fragile external dynamics [ Carlson and Doyle, 2002| . Similar robust-yet-fragile 



features have also been found in the of structure of diverse networks Peong et al., 2000| [Jeong et al., 2001 
[Albert etal, 20001 . 



Table 5: Comparison of synchronous and asynchronous algorithms. 





Synchronous 


Totally Asynchronous 


Random Order 


Two-timescale 


Assume 


Nodes are updated 


The time between 


Each node is 


In each time 




at multiples 


updates is 


updated at a randomly 


interval proteins are 




of the unit 


perturbed in a 


selected point of the 


updated first. 




time interval. 


range ite. 


unit time interval. 


then mRNAs. 


Update 


= k 




T'' = k-l + ^(t)^ 
cj)^ - node permutation 


= A; - 1 + 1^4, 


Pros 


Correctly identifies 


Allows for unlimited 


Does not depend 


Allows separation 




all steady states. 


variability in 


on any perturbation 


of post-translational 




Can be solved 


process durations. 


parameter e. 


and pre-translational 




analytically. 






processes. 


Cons 


Dynamics is 


Can have unrealistically short transcription 


Only useful when 




unrealistic. 


and translation times. 




process durations 
can be separated. 


Results 


Prepattem errors 


Divergence from the wild type process is possible. 


Development is stable if 




can be corrected. 


Cause: imbalance between two transcription factors. 


PTC is prepatterned. 



All our algorithms concur in suggesting that the divergence from wild type can be attributed to an im- 
balance between the two opposing Cubitus Interruptus transcription factors (CIA, CIR) in the posterior half 
of the parasegment. Thus the complementary regulation and pattern of these opposing transcription factors 
(Aza-Blanc and Kornberg 1999) is a vital requirement for the correct functioning of the segment polarity 
gene network. The totally asynchronous algorithm predicts that perturbations to the post-translational mod- 
ification of Cubitus Interruptus can have effects as severe as mutations: a transient overexpression of CIR 
leads to the pattern with no segmentation, while transient expression of CIA and not CIR leads to the broad 
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striped pattern. With the two timescale algorithm we find that the condition for the divergence from the wild 
type pattern is that, in the third cell of the parasegment, the post-translational modification of CI precedes 
the synthesis of the Patched protein. The biological realization of this condition appears unhkely, since 
PTC is documented as being ubiquitously expressed during cellularization (stage 5) ^Taylor et al., 1993J , 
while the post-translational modification of CI requires SMO that is only weakly expressed until stage 
8 lAlcedo et al., 20 00 1. Our model predicts that if for any reason the PTC protein is absent in the period 
when the pair-rule proteins decay and the regulation between the segment polarity genes starts, the wild type 
expression pattern is unreachable." 

Our methods combine the benefits of discrete-state models with a continuum in timescales. In the 
absence of quantitative information, we considered every possible timescale or update order, but as the two- 
timescale model demonstrates, existing information can be easily incorporated. We were able to describe 
the system in a rigorous mathematical way, to identify the relatively few types of behavior possible in the 
system (the attractors in state space) and to theoretically prove the convergence toward these states. Our 
results underscore that predictive mathematical modeling is possible despite the scarcity of quantitative 
information on gene regulatory processes. 
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A Additional Proofs 

Proof of Fact \I] For T = 0, the statement follows directly from the model's equations and by using the 
assumptions A1-A2 repeatedly: 

WGl = 1, en\ = 1, ENf = 1, 

as well as 

d?'^ = 0, C/°'^'^ = 0, CIAI'^ = 0, CIRI'^ = 0, PrC7°'^ = 0. 
We have (using ^) 

ptc\ 
PTC{ 
CIR{ 

so that we conclude 

ptc°'^'^ = 0, PTCl'^ = 0, C7/i?2'3 = 0. 
We next prove the Fact by induction. First note that, for any t > 0: 

wg\ = l WGf^ = 1 ^ en*i+^ = 1 ^ EN{+^ = 1 



= CIA\ and not CIR\ and not EN{ 
= ptc\-^ or [PTCl" and not i^i^al, 
= Clf" and PTCl" and not HHI'^ and not hhi-\ 
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and this implies 

d*+2 = ^ = 0. (18) 

Now assume that the Fact holds for some T > 1 and that 

wgl = l, 0<t<T. 
By the induction hypothesis, we know that 

Cll = 0, 3<t<T + 2, and CIR\ = 0, < t < T + 2, 

By (CHl, wgj = 1 imphes C/^+^ = 0, and this together with C/^^+^ = also guarantees that CIR^^^ = 
0, as we wanted to show. □ 
Proof of Fact^ To prove part (a), assume that wgl = for all < t < T, with T > 1. Since erig = 0, 
then EN2' = and hh^ = 0, HH2' = 0. For t > 3 apply Fact|5]to obtain the desired value for /1/12 and 
HH2. Note that ptcl = 1 implies PTCl = 1 and, together with HHI'"^ = 0, also PTCf = 1; then the 
value of PTC^ follows from Fact|5] For T > 2, and using ^ and ((SJl), we have 

CIAl = not PTCl" ov H Hi'' ovhh*^^ 
CIR\ = PTCl" and not HH'^' and not hh^^^ 

so, the values for /1/12, HH2 and PTC3, indeed imply that CIA\ = and CIR\ = 1, for 2 < i < T + 3. 

To prove part (b), note that if ci^ = 0, then also cig ^ = and hence CIA\ = 0. But now CIA\ = 
together with wg^, = immediately imply that wg\ = 0, and therefore, the results in part (a) are valid for 
all T > 0. □ 



Fact 4. 



/i/i^+2 = = and notCI/^^+^ for all T > 0. 



In particular, if C/i?^ = for all t > 3, then hh^^"^ = HHJ+^ = wg^, for all T > 1. 
Proof. Given any t > it is easy to see that 

WGl+' = wgl 
en*2+^ = WG\+^ or = H^Gg+S 

= £;iV*+2 and not C/i^*2"*"^ 

where the equation for en^2^^ follows from Q. □ 

Fact 5. PTCj = 1 and PrCj+^ = 0, for some T > 0, only if wgl = 1 for some t G {T - 3, T - 2} 
(chosen according to the permutation 

Proof. To see this, simply notice that 

PTCl""^ = ptcl or [PTCj and not i^i?*'' and not HH^^] 
= ptcl or [PTCj and not iJi?*''] 
= ptcl or [PrCj( and not uic/g""^ or C/Pg""^] 

because from ^ HHj = 0, and by Fact|4] Note that ta G {T, T + l}, depending on the permutation (f^'^^. 
So, for PTC3 to vanish it is necessary that both ptcj = and wg'^^~^ = 1. □ 
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B Attractiveness of the wild type pattern 



Assuming that the trajectory is not of the form ( fTTl . the only accessible steady state is the wild type. In 
this case, to establish convergence of the trajectory, it is enough to show that each node attains a constant 
value after a finite number of iterates. And in fact, from Propositions 14.11 14.41 and Corollaries 14.21 14.31 all 
the nodes become fixed after at most t iterates, as indicated: 

wgi,2 = 0, WGi^2 = 0, t>0, (H 

wg4 = 1, = 1, t>l, Proposition HTTI 

wgs = 0, WG3 = 0, t>l, Proposition E31 

eni = 1, ENi = 1, t>2, CorollarvIO 

en2 = 0, EN2 = 0, t>2, CorollarvIO 

en3,4 = 0, SA^3,4 = 0, i > 0, ^ 

hhi = 1, HHi = 1, t>3, CorollarvIO 

hh2 = 0, HH2 = 0, t>3, CorollarylO 

/l/l3,4 = 0, HHs,4 = 0, t>0, © 







cil - 


= 0, 


Ch = 


0, 


t>3, Corollary O 






ci2 -- 


= 1, 


CI2 = 


1, 


t>3, Corollary 10 






cisA = 


1, C/3,4 = 


1, 


t>2, (0, © 


ptci 


= 0, 


PTCi 


= 0, 


t > 0, 




Corollary 


CIAi 


= 0, 


CIRi 


= 0, 


t > 4, 




Corollary O 


CIAs 


= 0, 


CIR3 


= 1, 


t > 2, 




FactEl 


ptcs 


= 0, 


PTCs 


= 1, 


, t > 1, 




FactE 


CIA2 


= 1, 


CIR2 


= 0, 


t > 4, 




CI2 = 1 


ptC2 


= 1, 


PTC2 


= 1, 


t > 5, 




CIA2 = 1, CIR2 = 0, 


CIAi 


= 1, 


CIRi 


= 0, 


t > 3, 




wg4 = 1 


ptCi 


= 1, 


PTCi 


= 1, 


t > 5, 




CIAi = 1, CIRi = 0. 



This is indeed a complete characterization of the wild type steady state. 

Proof of Proposition Let P, A, C and R (< L) denote the positions of PTC3, CI A3, CI3 and 
CIR3, respectively. Then, from (flTl it is easy to see that 

PG {4,5,6,. ..,L}, ^G {2,3,4,. C G {1, 2, 3, . . . , A - 1}, R {I, . . . , P - l}\ {P, A,C}. 

To derive a formula for M, we note that, for each pair of values P, A, the number of possible combinations 
of C and R is: 

A — R P: {A-l){P-l-A) 



— R A P: {A-l){A-2), 
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respectively for sequences of the form CARP (top), or CRAP and RCAP (bottom). Therefore, summing 
over all posible A and P: 

: ^^'p-l)(P-l-^) + (^-l)(^-2)] 



Mr. 



p=4 A=2 
L P-l 



5^5^(^-l)(P-3) 

p=4 A=2 

}-J2{P-3){P-2){P-l) 

i£'i(i + i)(i + 2). 



Now, for L = 4, 



(L-4)! ^ 1 3, 0! ^ 1 1 ^ 1 
L! 2 ' 4! 2 4 8' 



Assume now that the equality is true for L — 1: 

M^_i(L-5)! 1 



Then 



(L-1)! 



= M^_-, + -(L-3)(L-2)(L-l) 

= i(L - 1)(L - 2)(L - 3)(L - 4) + {h - 3){L - 2){L - 1) 



= (L-l)(L-2)(L-3) 

= (L-l)(L-2)(L-3) 
just as we wanted to show. 



L 



1 L! 



8 (L-4)!' 



□ 



C State aggregation in the Markov chain 

The tables below show the complete transition probabilities pij (when not indicated, the transition probabil- 
ities are equal to 1). The states numbered 3 and 44 denote, respectively, the wild type and the mutant state. 
The initial condition was numbered 48. The first table shows the complete transition probabihties at step 1, 
from the wild type initial condition. 

Thus the probability shown in the diagram for the transition from the initial state to the aggregated state 
was obtained by adding p48,i + P48 6 = 0.1667 + 0.0417 = 0.2084. A more complex aggregation 



1 6 



formula was used for the transition 



10 12 



16 22 



1/ ^ 1/ 

2(^10,16+^10,22) + 2(^12,16 +P12,22) 



= - (0.2083 + 0.0417 + 0.3 + 0.0333) = 0.29165 « 0.29. 
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The normalization by 1/2 is justified by the fact that the transition from 9 and 30 to either 10 or 12 is the 
same. 



From initial wild type state to: 


1 


2 


3 


4 


5 


6 


0.1667 


0.357 


0.1667 


0.125 


0.125 


0.0417 



From 9 to: 


10 


11 


72 


13 


0.25 


0.25 


0.25 


0.25 





7 


16 


18 


20 


27 


22 


25 


26 


From 70 to: 


0.2083 


0.2083 


0.125 


0.0417 


0.125 


0.0417 


0.125 


0.125 


From 72 to: 


0.1167 


0.3 


0.1167 


0.05 


0.1333 


0.0333 


0.2 


0.05 





14 


15 


17 


19 


23 


24 


27 


28 


From ii to: 


0.2083 


0.125 


0.125 


0.2083 


0.125 


0.0417 


0.125 


0.0417 


From 13 to: 


0.1167 


0.1167 


0.1333 


0.3 


0.2 


0.0333 


0.05 


0.05 



From 29 or 31 to: 


32 


33 


34 


J5 


0.25 


0.25 


0.25 


0.25 



From 30 to: 


iO 


ii 


12 


ii 


0.25 


0.25 


0.25 


0.25 



From 36 to: 


5i 


37 


5S 


59 


40 


4i 


42 


43 


0.2083 


0.125 


0.125 


0.2083 


0.125 


0.125 


0.0417 


0.0417 



From 37 or 38 or 42 to: 


52 


33 


54 


55 


0.25 


0.25 


0.25 


0.25 



From 59 or 41 to: 




46 


0.5 


0.5 



From 40 or 45 to: 


44 


45 


46 


47 


0.25 


0.25 


0.25 


0.25 



References 



[Albert et al, 2000] Albert, R., Jeong, H. & Barabasi, A.-L., 2000. Error and attack tolerance in complex 
networks. Nature, 406, 378-382. 

[Albert and Othmer, 2003] Albert, R. & Othmer, H. G., 2003. The topology of the regulatory interactions 
predicts the expression pattern of the Drosophila segment polarity genes. J. Theor. Biol. 223, 1-18. 



26 



[Alcedo et al., 2000] Alcedo, J., Zou, Y. & Noll, M., 2000. Posttranscriptional regulation of smoothened is 
part of a self-correcting mechanism in the hedgehog signaling System. Molecular Cell 6, 457-465. 

[Alon et al., 1999] Alon, U., Surette, M., Barkai, N. and Leibler, S., 1999. Robustness in Bacterial Chemo- 
taxis. Nature 397, 168-171. 

[Aza-Blanc & Komberg 1999] Aza-Blanc, P. & Komberg, T. B. (1999) Ci, a complex transducer of the 
Hedgehog signal. Trends in Genetics 15, 458-462. 

[Bejsovec and Wieschaus, 1993] Bejsovec, A. & Wieschaus, E., 1993. Segment polarity gene interactions 
modulate epidermal patterning in Drosophila embryos. Development 119, 501-517. 

[Bemot et al, 2004] Bernot, G., Comet, J. -P., Richard, A. & Guespin, J., 2004. Application of formal meth- 
ods to biological regulatory networks: extending Thomas' asynchronous logical approach with tem- 
poral logic. J. Theor. Biol. 229, 339-347. 

[Bertsekas and Tsitsiklis, 1989] Bertsekas, D.P. & Tsitsiklis, J.N., 1989. Parallel and Distributed Computa- 
tion, Numerical Methods. Prentice Hall, New Jersey. 

[Bodnar, 1997] Bodnar, J. W., 1997. Programming the Drosophila Embryo. J. Theor. Biol. 188, 391-445. 

[Cadigan et al., 1994] Cadigan, K. M., Grossniklaus, U. & Gehring, W. J., 1994. Localized expression of 
sloppy paired protein maintains the polarity of Drosophila parasegments. Genes & Dev. 8, 899-913. 

[Cadigan & Nusse 1997] Cadigan, K. M., & Nusse, R. (1997) Wnt signaling: a common theme in animal 
development. Genes Dev. 11, 3286-3305. 

[Carlson and Doyle, 2002] Carlson, J. M. & Doyle, J., 2002. Complexity and robustness. Proc. Nat. Acad. 
Sci. 99, 2538-2545. 

[Conant and Wagner, 2004] Conant, G. C.& Wagner A., 2004. Duplicate genes and robustness to transient 
gene knock-downs in Caenorhabditis elegans. Proc. R. Soc. Lond. B Biol. Sci. 271, 89-96. 

[Davidson et al., 2002] Davidson, E. H., et al., 2002. A genomic regulatory network for development. Sci- 
ence 295, 1669-1678. 

[de Jong et al., 2004] de Jong, H. et al. , 2004. Qualitative simulation of genetic regulatory networks using 
piecewise-linear models. Bull. Math. Bio. 66, 301-340. 

[DiNardo et al., 1988] DiNardo, S., Sher, E., Heemskerk-Jongens, J., Kassis, J. A. & O'Farrell, P. H., 1988. 
Two-tiered regulation of spatially patterned engrailed gene expression during Drosophila embryoge- 
nesis. Nature 332, 45-53. 

[Eaton & Komberg 1990] Eaton, S. & Komberg, T. B. (1990) Repression of ci-D in posterior compartments 
of Drosophila by engrailed. Genes. Dev. 4, 1068-1077. 

[Eldar et al., 2002] Eldar, A., Dorfman, R., Weiss, D., Ashe, H., Shilo, B.-Z. & Barkai, N., 2002. Robust- 
ness of the BMP morphogen gradient in Drosophila embryonic patterning. Nature 419, 304-308. 

[Feller, 1970] Feller, W., 1970. An Introduction to Probability Theory and its Applications (Third edition, 
revised). John Wiley & Sons, New York, 1970. 



27 



[Gallet et al., 2000] Gallet, A., Angelats, C, Kerridge, S. & Therond, P. P., 2000. Cubitus intemiptus- 
independent transduction of the Hedgehog signal in Drosophila, Development 127, 5509-5522. 

[Ghysen and Thomas, 2003] Ghysen, A. & Thomas, R., 2003. The formation of sense organs in Drosophila: 
A logical approach. BioEssays 25, 802-807. 

[Glass and Kauffman, 1973] Glass, L. & Kauffman, S. A., 1973. The logical analysis of continuous, non- 
linear biochemical control networks. J. Theor. Biol. 39, 103-129. 

[Grossniklaus et al., 1992] Grossniklaus, U., Pearson, R. K. & Gehring, W. J. (1992) The Drosophila sloppy 
paired locus encodes two proteins involved in segmentation that show homology with mammalian 
transcription iaciors.Genes Dev. 6, 1030-1051. 

[Gursky et al., 2001] Gursky, V. V., Reinitz, J. & Samsonov, A. M., 2001. How gap genes make their do- 
mains: An analytical study based on data driven approximations. Chaos 11, 132-141. 

[Hidalgo and Ingham, 1990] Hidalgo, A. & Ingham, P., 1990. Cell Patterning in the Drosophila segment: 
spatial regulation of the segment polarity gene patched. Development 110, 291-301. 

[Hooper and Scott, 1992] Hooper, J. E. & Scott, M. P., 1992. The Molecular Genetic Basis of Positional 
Information in Insect Segments. In: Early Embryonic Development of Animals (ed. Hennig, W.) 
1-49, Springer, Berlin. 

[Ingham 1998] Ingham, R W. (1998) Transducing hedgehog: the story so far, EMBO J. 17, 3505-3511. 

[Ingham & McMahon 2001] Ingham, P. W. & McMahon, A. P. (2001) Hedgehog signahng in animal de- 
velopment: paradigms and principles. Genes Dev. 15, 3059-3087. 

[Ingham et al., 1991] Ingham, P.W., Taylor, A. M. & Nakano, Y., 1991. Role of the Drosophila patched 
gene in positional signaling. Nature 353, 184-187. 

[Jeong et al, 2000] Jeong, H., Tombor, B., Albert, R., Oltvai, Z.N. & Barabasi, A.-L., 2000. The large-scale 
organization of metabolic networks. Nature 407, 651-654. 

[Jeong et al., 2001] Jeong, H., Mason, S., Barabasi, A.-L. & Oltvai, Z. N., 2001. Lethality and centrality in 
protein networks. Nature 411, 41-42. 

[Kauffman, 1993] Kauffman, S. A., 1993. The origins of Order. Oxford University Press, New York. 

[Kauffman et al., 2003] Kauffman, S., Peterson, C, Samuelsson, B., Troein, C, 2003. Random Boolean 
network models and the yeast transcriptional network. Proc. Natl. Acad. Sci. USA. 100 14796-9. 

[Martinez-Arias et al., 1988] Martinez-Arias, A., Baker, N. & Ingham, P. W., 1988. Role of segment polar- 
ity genes in the definition and maintenance of cell states in the Drosophila embryo. Development 
103, 157-170. 

[Mendoza et al., 1999] Mendoza, L., Thieffry, D. & Alvarez-Buylla, E. R., 1999. Genetic control of flower 
morphogenesis in Arabidopsis thaliana: a logical analysis. Bioinformatics 15, 593-606. 

[Ohlmeyer & Kalderon 1998] Ohlmeyer, J. T. & Kalderon, D. (1998) hedgehog stimulates maturation of 
Cubitus interruptus into a labile transcriptional activator. Nature 396, 749-753. 



28 



[Pfeiffer & Vincent 1999] Pfeiffer, S. & Vincent, J.-P. (1999) Signaling at a distance:Transport of Wingless 
in the embryonic epidermis of Drosophila. Cell & Dev. Biol. 10, 303-309. 

[Rao et al., 2002] Rao, C.V, Wolf, D.M. & Arkin, A. P., 2002. Control, exploitation and tolerance of intra- 
cellular noise. Nature 420, 231-237. 

[Reinitz and Sharp, 1995] Reinitz, J. & Sharp, D. H., 1995. Mechanism of eve stripe formation. Mecha- 
nisms of Development 49, 133-158. 

[Sanchez and Thieffry, 2001] Sanchez, L., & Thieffry, D., 2001. A logical analysis of the Drosophila gap- 
gene system. J. Theor. Biol. 211, 115-141. 

[Schwartz et al., 1995] Schwartz, C, Locke, J., Nishida, C. & Kornberg, T. B., 1995. Analysis of cubitus 
interruptus regulation in Drosophila embryos and imaginal disks. Development 121, 1625-1635. 

[Tabata et al, 1992] Tabata, T., Eaton, S. & Kornberg, T. B., 1992. The Drosophila hedgehog gene is 
expressed specifically in posterior compartment cells and is a target of engrailed regulation. Genes & 
Dev. 6, 2635-2645. 

[Taylor et al, 1993] Taylor, A. M., Nakano, Y, Mohler, J. & Ingham, P. W., 1993. Contrasting distributions 
of patched and hedgehog proteins in the Drosophila embryo. Mechanisms of Development 42, 89-96. 

[Thomas, 1973] Thomas, R., 1973. Boolean formalization of genetic control circuits. J. Theor. Biol. 42, 
563-585. 

[von Dassow et al., 2000] von Dassow, G., Meir., E., Munro, E. M., & Odell, G. M., 2000. The segment 
polarity network is a robust developmental module. Nature 406, 188-192. 

[Wolpert et al., 1998] Wolpert, L., Beddington, R., Brockes, J., Jessell, T., Lawrence, P., & Meyerowitz, E., 
1998. Principles of Development Current Biology Ltd., London. 

[Yuh et al, 2001] Yuh, C. H., Bolouri, H., Bower, J. M. and Davidson, E. H., 2001. A logical model of cis- 
regulatory control in a eukaryotic system. In: Computational Modeling of Genetic and Biochemical 
Networks (eds. Bower, J. M. and Bolouri, H.), 73-100, MIT Press, Cambridge, MA. 



29 



